to remove/mbnmadoseMODEL.R

In htx-r/doseresNMA:

doseresNMA <- function(){
			# Begin Model Code

    s.beta.2[ref.index] <- 0

    s.beta.1[ref.index] <- 0

    for(i in 1:ns){ # Run through all NS trials

      delta[i,1] <- 0

      #DR[i,1] <- 0 # Dose-response model is 0 for baseline arms
      mu[i] ~ dnorm(0,0.001)

      for (k in 1:na[i]){ # Run through all arms within a study

        rhat[i,k] <- p[i,k] * n[i,k]
        resdev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (n[i,k] - r[i,k]) * (log(n[i,k] - r[i,k]) - log(n[i,k] - rhat[i,k])))
        r[i,k] ~ dbin(p[i,k], n[i,k])
        logit(p[i,k]) <- theta[i,k]
        theta[i,k] <- mu[i] + delta[i,k]
      }

      resstudydev[i] <- sum(resdev[i, 1:na[i]])

      for(k in 2:na[i]){ # Treatment effects

        delta[i,k] <- DR[i,k]

        DR[i,k] <- ((s.beta.1[t[i,k]] * dose1[i,k]) + (s.beta.2[t[i,k]] * dose2[i,k])) - ((s.beta.1[t[i,1]] * dose1[i,1]) + (s.beta.2[t[i,1]] * dose2[i,1]))
      }
    }

    for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)){ # Priors on relative treatment effects
      mult[1:2,k] ~ dmnorm(d.prior[], inv.R[1:2, 1:2])
      d.2[k] <- mult[2,k]

      s.beta.2[k] <- d.2[k]

      d.1[k] <- mult[1,k]

      s.beta.1[k] <- d.1[k]
    }

    totresdev <- sum(resstudydev[])
    Omega[1,1] <- 1
    Omega[2,2] <- 1

    for (r in 1:2) {
      d.prior[r] <- 0
    }

    inv.R ~ dwish(Omega[,], 2)

    for (r in 1:(2-1)) {  # Covariance matrix upper/lower triangles
      for (c in (r+1):2) {
        Omega[r,c] <- 0   # Lower triangle
        Omega[c,r] <- 0   # Upper triangle
      }
    }
}



mydoseresNMA <- function(){
  b1[ref.index] <- 0
  b2[ref.index] <- 0
    for (i in 1:ns) {
      delta[i, 1] <- 0
      #d[i, 1] <- 0.00000E+00
      u[i] ~ dnorm(0, 0.001)
      for (k in 1:na[i]) {
        r[i, k] ~ dbin(p[i, k], n[i, k])
        logit(p[i, k]) <- theta[i,k]
        theta[i,k] <- u[i] + delta[i, k]
        rhat[i, k] <- p[i, k] * n[i, k]
        dev[i, k] <- 2 * (r[i, k] * (log(r[i, k]) - log(rhat[i,k])) + (n[i,k] - r[i, k]) * (log(n[i, k] - r[i,k]) - log(n[i, k] - rhat[i, k])))
      }
      resdev[i] <- sum(dev[i, 1:na[i]])
      for (k in 2:na[i]) {
        d[i, k] <- (b1[t[i, k]] * dose1[i, k]) - (b1[t[i, 1]] * dose1[i, 1]) + (b2[t[i, k]] *dose2[i, k]) - (b2[t[i, 1]] * dose2[i,1])
        delta[i, k] <- d[i, k]
      }
    }

    for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)) {
      # b[1:2,k]~dmnorm(b.prior[],inv.R[1:2,1:2])
      # # b1[k] <- b[1,k]
      # # b2[k] <- b[2,k]
      # d1[k] <- b[1,k]
      # d2[k] <- b[2,k]
      #
      # b1[k] <- d1[k]
      # b2[k] <- d2[k]

      b1[k] ~ dnorm(0.00000E+00, 0.001)
      b2[k] ~ dnorm(0.00000E+00, 0.001)
    }

    # b.prior[1] <- 0
    # b.prior[2] <- 0
  for (r in 1:2) {
    b.prior[r] <- 0
  }


    inv.R ~ dwish(Omega[,], 2)
    Omega[1,1] <- 1
    Omega[2,2] <- 1

    # Omega[1,2] <- 0
    # Omega[2,1] <- 0

    for (r in 1:(2-1)) {  # Covariance matrix upper/lower triangles
      for (c in (r+1):2) {
        Omega[r,c] <- 0   # Lower triangle
        Omega[c,r] <- 0   # Upper triangle
      }
    }



}


dosresnetmeta.jagsdata <- makeJAGSdata(data=antidep,metareg=F,class.effect =F,ref.lab='placebo',refclass.lab = 'placebo')
dosresnetmeta.jagsdata$rmin <- dosresnetmeta.jagsdata$r[,1]
dosresnetmeta.jagsdata$nmin <- dosresnetmeta.jagsdata$n[,1]

# jags model
#jagsmodelDRnetmeta<- makeJAGSmodel(metareg=F,beta.effects='random',consistency=T,class.effect=F,doselink='spline')
# run jags model
doseresNMAfit <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('s.beta.1','s.beta.2'),model.file = doseresNMA,
                                          n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)

histDRparam2(doseresNMAfit)


# 7. Histogram of DR parameters
histDRparam2 <- function(x=x,metareg=F,class.effect=F){ # ,  drug.lab=drug.lab
  require(tidyr)

  param.lab <- c('s.beta.1.','s.beta.2.')

  if(metareg){
    param.lab0 <- 'g.'
    param.lab <- append(param.lab,param.lab0)
  }
  sims <- x$BUGSoutput$sims.array
  sims.mat <- apply(sims,3L,c) # merge the 3 chains (columns into one)
  sims.param <- sims.mat%>%data.frame()%>%select(starts_with(param.lab))
  plotdata <- data.frame(sim=unlist(c(sims.param)),param=rep(colnames(sims.param),each=nrow(sims.param)))

  g <- ggplot(plotdata, aes(x=sim)) +
    geom_histogram(alpha=0.6, binwidth = 0.001,color="darkblue", fill="lightblue")
  g <- g + ggplot2::facet_wrap(~plotdata$param,scales = 'free')
  g <- g + ggplot2::labs(y="Frequency", x="iterations")
  g+ theme(axis.line = element_line(colour = "black"),
           panel.grid.major = element_blank(),
           panel.grid.minor = element_blank())+theme_bw()
}



jagsmodelDRnetmeta<- makeJAGSmodel2(metareg=F,beta.effects='fixed',consistency=T,class.effect=F,doselink='spline')
# run jags model
dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('b1','b2'),model.file = jagsmodelDRnetmeta,
                                          n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
dosresnetmeta.runjagsRCS$model

mydoseresNMAfit <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save =  c('s.beta.1','s.beta.2'),model.file = mydoseresNMA,
                                          n.chains=2,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 5)
mydoseresNMAfit$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),]

res12 <- round(cbind(doseresNMAfit$BUGSoutput$mean$s.beta.1,dosresnetmeta.runjagsRCS$BUGSoutput$mean$b1),4)
dif <- unlist(sapply(1:nrow(res12), function(i) res12[i,1]-res12[i,2]))
cbind(res12,dif=dif)

levels(factor(antidep$drug))
antidep[antidep$drug=='reboxetine',]

round(cbind(doseresNMAfit$BUGSoutput$mean$s.beta.2,dosresnetmeta.runjagsRCS$BUGSoutput$mean$b2,mydoseresNMAfit$BUGSoutput$mean$b2),4)
round(rbind(doseresNMAfit$BUGSoutput$summary[c('s.beta.1[17]','s.beta.2[17]'),],dosresnetmeta.runjagsRCS$BUGSoutput$summary[c('b1[17]','b2[17]'),],mydoseresNMAfit$BUGSoutput$summary[c('b1[17]','b2[17]'),]),4)

# dose ditsribution per drug
antidep$tindex <- as.factor(as.numeric(factor(antidep$drug)))
ggplot(data = antidep, aes(x=tindex, y=dose)) +
  facet_wrap( ~ drug, scales="free")+
  geom_dotplot(binaxis='y', stackdir='center', dotsize=3,col='orange')